function [] = runCounterfactuals(runId, cnfactType)

rng(1); % for replicability

%%% Add paths
addpath('../general_funcs');
addpath('../model_funcs');
addpath('../model_funcs/dim_specific_funcs/dim2');
addpath('../model_funcs/dim_specific_funcs/dim3');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');

%%% Get input / output paths
config = getConfig_jointProcess(runId);

% Get config_consumerAdoptions and config_providerAdoptions: necessary input for make_Xdynamic_xxx functions
indivConfigs.consumerAdoptions = getConfig_consumerAdoptions(config.runId_consumerAdoptions);
indivConfigs.providerAdoptions  = getConfig_providerAdoptions(config.runId_providerAdoptions);
indivConfigs.providerExits      = getConfig_providerExits(config.runId_providerExits);

indivRunFolders.providerAdoptions  = config.runFolder_providerAdoptions;
indivRunFolders.providerExits      = config.runFolder_providerExits;
indivRunFolders.consumerAdoptions = config.runFolder_consumerAdoptions;

cnfactName = sprintf('cnfact%d', cnfactType);

outputFolder         = createFolderIfNotExist(sprintf('%s/cnfacts/%s', config.runFolder, cnfactName));
baseConditionsFolder = sprintf('%s/dynamicPredictions/periodStartConditions', config.runFolder);
basePredMatFile2     = sprintf('%s/../pred_vals.mat', baseConditionsFolder);

outputMatFile = sprintf('%s/cnfact_vals.mat', outputFolder);

%%% Load data from the 3 submodels
data = load_data_fullmodel(indivRunFolders);

%%% Load Xdynamic_staticInputs
[~, Xdynamic_staticInputs] = load_base_initial_conditions_and_Xdynamic_staticInputs(indivRunFolders);

%%%  Load parameter estimates from the 3 submodels
[~,paramsParts] = load_param_estimates(data, indivRunFolders);

%%% Aggregate the data from spell level to day level to get day-level covariates
data.consumerAdoptions = aggregateCtsTimeData_consumerAdoptions(data.consumerAdoptions);
data.providerAdoptions  = aggregateCtsTimeData_providerAdoptions(data.providerAdoptions);
data.providerExits      = aggregateCtsTimeData_providerExits(data.providerExits);

%%% Read dimensions
T  = data.consumerAdoptions.dims.dims1{4};
K1 = data.consumerAdoptions.dims.dims1{5};
K2 = data.consumerAdoptions.dims.dims1{6};
K  = data.providerAdoptions.dims.dims1{3};

%%% Combine dims into object mydims
mydims.T  = T;
mydims.K  = K;
mydims.K1 = K1;
mydims.K2 = K2;

% Get labels
periodLabels = data.periodLabels;
geoLabels    = data.geoLabels;

% Get geo fields useful to define counterfactual scenarios
population        = data.population; % K x 1
surfaces          = data.surfaces; % K x 1
incomingCommuters = data.incomingCommuters; % K x 1
touristBeds       = data.touristBeds; % K x 1

%%%%%%%%%%%%%%%%%% DEFINE COUNTERFACTUAL SCENARIO %%%%%%%%%%%%%%%%%%
if cnfactType == 1
	cnfactLabel = 'Only consumers seeded'; tt_intervention = 1;
end
if cnfactType == 2
	cnfactLabel = 'Only providers seeded'; tt_intervention = 1;
end

if cnfactType == 3
	cnfactLabel = 'Seeding providers in top 5 cities'; tt_intervention = 1;
end
if cnfactType == 4
	cnfactLabel = 'Seeding providers proportional to tourist lodging capacity'; tt_intervention = 1;
end
if cnfactType == 5
	cnfactLabel = 'Seeding providers proportional to population'; tt_intervention = 1;
end

if cnfactType == 6
	cnfactLabel = 'Remove Paris subway campaign';
	myVarName = 'photoContest_IDF';
	
	% Set dummy equal to zero in consumerAdoptions covariates
	disp(sprintf('Setting the following X equal to zero in consumerAdoptions: %s', myVarName));
	idx = find(strcmp(Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.Xnames, myVarName));
	if isempty(idx); error(sprintf('Could not find %s in X_tk1 names for consumerAdoptions', myVarName)); end;
	
	% Make tt_intervention
	tmp = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X(:,idx);
	tmp = reshape(tmp, [T K1]);
	tt_intervention = min(find(any(tmp > 0, 2)));
	disp(sprintf('tt_intervention: %d', tt_intervention));
	clear tmp;
	% Set dummy equal to zero
	Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X(:,idx) = 0;
	
	% Set dummy equal to zero in providerAdoptions covariates
	disp(sprintf('Setting the following X equal to zero in providerAdoptions: %s', myVarName));
	idx = find(strcmp(Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.Xnames, myVarName));
	if isempty(idx); error(sprintf('Could not find %s in X_tk names for providerAdoptions', myVarName)); end;
	Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X(:,idx) = 0;
end


%%% Compute the part of V that does not depend on installedBases
XbetaFixed = calc_Xbeta_fixed(data, paramsParts);


%%% Preprocess Xdynamic_staticInputs inputs to compute the parts of V that depend on installedBases
Xdynamic_staticInputs = permute_Xparts_static_X_X_FEs(Xdynamic_staticInputs, mydims);


%%% Subset XbetaFixed and Xdynamic_staticInputs: retain periods from tt_intervention to T (beforehand: everything is the same as in the base case)
XbetaFixed.providerAdoptions  = XbetaFixed.providerAdoptions(:,tt_intervention:end); % K x T_residual
XbetaFixed.providerExits      = XbetaFixed.providerExits(:,tt_intervention:end); % K x T_residual
XbetaFixed.consumerAdoptions = XbetaFixed.consumerAdoptions(:,:,tt_intervention:end); % K x T_residual

Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X     = Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X(:,:,tt_intervention:end);
Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X_FEs = Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X_FEs(:,:,tt_intervention:end);

Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X     = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X(:,:,tt_intervention:end);
Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X_FEs = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{2}.X_FEs(:,:,tt_intervention:end);

Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X     = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X(:,:,tt_intervention:end);
Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X_FEs = Xdynamic_staticInputs.consumerAdoptions.Xparts_static{3}.X_FEs(:,:,tt_intervention:end);

% Get (final) dimensions
T = size(Xdynamic_staticInputs.providerAdoptions.Xparts_static{1}.X, 3);


%%%% Launch predictions in counterfactual scenario
NumPreds = 100;
cnfact_vals.consumerAdoptions   = zeros(T,K1,NumPreds); % it is Y_tk1 that is going to be stored (too much storage required for Y_tk1tk2)
cnfact_vals.providerAdoptions    = zeros(T,K,NumPreds);
cnfact_vals.providerExits        = zeros(T,K,NumPreds);
cnfact_vals.laggedProviders      = zeros(T,K,NumPreds);
cnfact_vals.laggedConsumersOrig = zeros(T,K,NumPreds);

tic
for pp = 1:NumPreds
	% From base case (dynamic predictions): load the initial conditions at start of period tt_intervention for sequence draw pp
	load(sprintf('%s/periodStartConditions_%d.mat', baseConditionsFolder, pp), 'periodStartConditions');
	initConds_pp.laggedProviders = periodStartConditions{tt_intervention}.laggedProviders;
	initConds_pp.laggedConsumersOrig = periodStartConditions{tt_intervention}.laggedConsumersOrig;
	initConds_pp.laggedConsumersDest = periodStartConditions{tt_intervention}.laggedConsumersDest;
	initConds_pp.popAtRisk = periodStartConditions{tt_intervention}.popAtRisk;
	initConds_pp.randomState = periodStartConditions{tt_intervention}.randomState;
	rng(initConds_pp.randomState); % Set random seed where it was in base case at the start of the same period tt.
	clear periodStartConditions;
	
	%%%%% Make intervention on installed base at tt_intervention %%%%%
	if cnfactType <= 5
		aggLaggedProviders = sum(initConds_pp.laggedProviders);
		aggLaggedConsumers = sum(initConds_pp.laggedConsumersOrig);
				
		if cnfactType == 1 % Only consumers seeded
			consumerOrigChanges = zeros(K,1);
			consumerDestChanges = zeros(K,1);
			providerChanges      = - initConds_pp.laggedProviders;
		end
		
		if cnfactType == 2 % Only providers seeded
			consumerOrigChanges = - initConds_pp.laggedConsumersOrig;
			consumerDestChanges = - initConds_pp.laggedConsumersDest;
			providerChanges      = zeros(K,1);
		end
		
		if cnfactType == 3 % Providers distributed in top 5 cities only, no initial consumers
			NumCities = 5;
			[~, cities_idxes] = mink(-population, NumCities);
			mycriterion = zeros(K,1); mycriterion(cities_idxes) = 1;
			newLaggedProviders      = distribute_proportionally_2_integers(aggLaggedProviders, mycriterion);  % K x 1
			newLaggedConsumersOrig = zeros(K,1);
			newLaggedConsumersDest = zeros(K,1);
		end
		if cnfactType == 4 % Providers distributed prop. to tourist lodging capacity, no initial consumers
			newLaggedProviders      = distribute_proportionally_2_integers(aggLaggedProviders, touristBeds);  % K x 1
			newLaggedConsumersOrig = zeros(K,1);
			newLaggedConsumersDest = zeros(K,1);
		end
		if cnfactType == 5 % Providers distributed prop. to population, no initial consumers
			newLaggedProviders      = distribute_proportionally_2_integers(aggLaggedProviders, population);  % K x 1
			newLaggedConsumersOrig = zeros(K,1);
			newLaggedConsumersDest = zeros(K,1);
		end
		
		if cnfactType >= 3 && cnfactType <= 5
			providerChanges      = newLaggedProviders      - initConds_pp.laggedProviders;
			consumerOrigChanges = newLaggedConsumersOrig - initConds_pp.laggedConsumersOrig;
			consumerDestChanges = newLaggedConsumersDest - initConds_pp.laggedConsumersDest;
		end
		
		% Update lagged base + pop at risk based on changes made
		initConds_pp.laggedProviders      = initConds_pp.laggedProviders      + providerChanges;
		initConds_pp.laggedConsumersOrig = initConds_pp.laggedConsumersOrig + consumerOrigChanges;
		initConds_pp.laggedConsumersDest = initConds_pp.laggedConsumersDest + consumerDestChanges;
		initConds_pp.popAtRisk.providerAdoptions  = initConds_pp.popAtRisk.providerAdoptions  - providerChanges;
		initConds_pp.popAtRisk.providerExits      = initConds_pp.popAtRisk.providerExits      + providerChanges;
		initConds_pp.popAtRisk.consumerAdoptions = initConds_pp.popAtRisk.consumerAdoptions - consumerOrigChanges;
	end
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	[Ycnfact_pp, periodStartConditions_pp] = generate_fullmodel_sequence_ctsTime(indivConfigs, initConds_pp, XbetaFixed, Xdynamic_staticInputs, paramsParts);
	Ycnfact_pp.consumerAdoptions = reshape(Ycnfact_pp.consumerAdoptions, [T*K1 K2]); % (T*K1) x K2
	Ycnfact_pp.consumerAdoptions = reshape(sum(Ycnfact_pp.consumerAdoptions, 2), [T K1]); % T x K1

	% Get installed bases from Ycnfact
	Ycnfact_pp = compute_laggedBases(Ycnfact_pp, initConds_pp, false);
	
	% Store everything
	cnfact_vals.consumerAdoptions(:,:,pp)   = Ycnfact_pp.consumerAdoptions;
	cnfact_vals.providerAdoptions(:,:,pp)    = Ycnfact_pp.providerAdoptions;
	cnfact_vals.providerExits(:,:,pp)        = Ycnfact_pp.providerExits;
	cnfact_vals.laggedProviders(:,:,pp)      = Ycnfact_pp.laggedProviders;
	cnfact_vals.laggedConsumersOrig(:,:,pp) = Ycnfact_pp.laggedConsumersOrig;
	
	% Display time
	displayRemainingTime(pp, NumPreds, 1);
end
toc


%%%%% Load base case time series for periods before intervention
load(basePredMatFile2, 'pred_vals');
prev_vals.consumerAdoptions   = pred_vals.consumerAdoptions(1:tt_intervention-1,:,:);   % Tbefore x K x NumPreds
prev_vals.providerAdoptions    = pred_vals.providerAdoptions(1:tt_intervention-1,:,:);    % Tbefore x K x NumPreds
prev_vals.providerExits        = pred_vals.providerExits(1:tt_intervention-1,:,:);        % Tbefore x K x NumPreds
prev_vals.laggedProviders      = pred_vals.laggedProviders(1:tt_intervention-1,:,:);      % Tbefore x K x NumPreds
prev_vals.laggedConsumersOrig = pred_vals.laggedConsumersOrig(1:tt_intervention-1,:,:); % Tbefore x K x NumPreds
clear pred_vals;


%%%%% Combine before- and and after-intervention periods
[Tafter,K,NumPreds] = size(cnfact_vals.consumerAdoptions);
cnfact_vals.consumerAdoptions   = reshape(cnfact_vals.consumerAdoptions,   [Tafter K*NumPreds]); % Tafter x (K*NumPreds)
cnfact_vals.providerAdoptions    = reshape(cnfact_vals.providerAdoptions,    [Tafter K*NumPreds]); % Tafter x (K*NumPreds)
cnfact_vals.providerExits        = reshape(cnfact_vals.providerExits,        [Tafter K*NumPreds]); % Tafter x (K*NumPreds)
cnfact_vals.laggedProviders      = reshape(cnfact_vals.laggedProviders,      [Tafter K*NumPreds]); % Tafter x (K*NumPreds)
cnfact_vals.laggedConsumersOrig = reshape(cnfact_vals.laggedConsumersOrig, [Tafter K*NumPreds]); % Tafter x (K*NumPreds)

[Tbefore,K,NumPreds] = size(prev_vals.consumerAdoptions);
prev_vals.consumerAdoptions   = reshape(prev_vals.consumerAdoptions,   [Tbefore K*NumPreds]); % Tbefore x (K*NumPreds)
prev_vals.providerAdoptions    = reshape(prev_vals.providerAdoptions,    [Tbefore K*NumPreds]); % Tbefore x (K*NumPreds)
prev_vals.providerExits        = reshape(prev_vals.providerExits,        [Tbefore K*NumPreds]); % Tbefore x (K*NumPreds)
prev_vals.laggedProviders      = reshape(prev_vals.laggedProviders,      [Tbefore K*NumPreds]); % Tbefore x (K*NumPreds)
prev_vals.laggedConsumersOrig = reshape(prev_vals.laggedConsumersOrig, [Tbefore K*NumPreds]); % Tbefore x (K*NumPreds)

cnfact_vals.consumerAdoptions   = [prev_vals.consumerAdoptions;   cnfact_vals.consumerAdoptions];   % T x (K*NumPreds)
cnfact_vals.providerAdoptions    = [prev_vals.providerAdoptions;    cnfact_vals.providerAdoptions];    % T x (K*NumPreds)
cnfact_vals.providerExits        = [prev_vals.providerExits;        cnfact_vals.providerExits];        % T x (K*NumPreds)
cnfact_vals.laggedProviders      = [prev_vals.laggedProviders;      cnfact_vals.laggedProviders];      % T x (K*NumPreds)
cnfact_vals.laggedConsumersOrig = [prev_vals.laggedConsumersOrig; cnfact_vals.laggedConsumersOrig]; % T x (K*NumPreds)

T = Tbefore+Tafter;
cnfact_vals.consumerAdoptions   = reshape(cnfact_vals.consumerAdoptions,   [T K NumPreds]); % T x K x NumPreds
cnfact_vals.providerAdoptions    = reshape(cnfact_vals.providerAdoptions,    [T K NumPreds]); % T x K x NumPreds
cnfact_vals.providerExits        = reshape(cnfact_vals.providerExits,        [T K NumPreds]); % T x K x NumPreds
cnfact_vals.laggedProviders      = reshape(cnfact_vals.laggedProviders,      [T K NumPreds]); % T x K x NumPreds
cnfact_vals.laggedConsumersOrig = reshape(cnfact_vals.laggedConsumersOrig, [T K NumPreds]); % T x K x NumPreds


%%%% Aggregate across locations
cnfact_aggVals = calc_Yagg(cnfact_vals);

%%% Save predictions to output file 1
save(outputMatFile, 'cnfact_vals', 'cnfact_aggVals', 'cnfactName', 'cnfactLabel', 'geoLabels', 'periodLabels', '-v7.3');

end